clear all
close all
filename = 515;
[Vdc, Tset, T_sample, I2,~,Vxx,~] = ldsm(getFilenameFromIndex(filename));
Vdc = Vdc+1.3e-4;

filename = 516;
[Vdc2, Tset2, T_sample2, I22,~,Vxx2,~] = ldsm(getFilenameFromIndex(filename));
Vdc2 = Vdc2+1.3e-4;

filename = 517;
[Vdc3, Tset3, T_sample3, I23,~,Vxx3,~] = ldsm(getFilenameFromIndex(filename));
Vdc3 = Vdc3+1.3e-4;

Vdc = cat(2, Vdc, Vdc2, Vdc3);
T_sample = cat(2, T_sample, T_sample2, T_sample3);
I2 = cat(2, I2, I22, I23);
Vxx = cat(2, Vxx, Vxx2, Vxx3);
% The above data are an example for Dynes formula fitting fot Vtg = -4.55, 
% corresponding to the main Fig. 3B. The rest of the data to reproduce other 
% fittings in the manuscript can be found in our data base.


%% Superconducting Gap Fitting Script
% This script analyzes tunneling spectroscopy data by fitting it to
% Dynes-like models (nodal or nodeless) with Gaussian broadening.
% Outputs include fitting results, parameter extraction, and visualizations.

% INPUT DATA:
%   - I2, Vxx: m-by-n matrices (m = bias voltage points, n = temperatures)
%   - T_sample: m-by-n matrix of temperature values
%   - Vdc:  m-by-n matrices of bias voltage array
%   - Functions required: 
%       DynesFit_nodal_gauss_fraction_Base.m
%       DynesFit_swave_gauss_fraction_Base.m

%% Data Preprocessing
N = 21;                                  % N = 21 corresponds to T = 1.4K as the normaliztion curve
dIdV = I2 ./ Vxx;                        % Differential conductance
dIdV_nor = dIdV ./ dIdV(:,N);            % Normalize using the tunneling curve at the N-th temperature
T_ave = mean(T_sample, 1);               % Average temperature at each point

lineStyles = linspecer(length(T_ave),'sequential'); % Color map for plots

%% Analysis Parameters
Vtglist = [-4.55];        % Gate voltage at which the tunneling data was taken
fitRange = 1;             % Number of temperature traces to fit
xRange = [-0.5 0.5];      % Fitting range in meV
isNodal = true;           % Fitting option: true = d-wave, false = s-wave

% Initial fitting parameters:
%   [delta0, gamma, slope_BG, offset_BG, sigma]
param0 = [0.1, 0.02, 0.1, 1, 0.1];
ParamLowerBound = [0, 1e-5, -0.5, 0.5, 0];
ParamUpperBound = [0.5, 1, 0.5, 1.5, 1];

options = optimset('Display','iter','PlotFcns','optimplotfval');

%% Fitting Loop
clear fitParam
for i = 1:length(Vtglist)                % Loop over gate voltages
    for j = 1:fitRange                   % Loop over temperature traces
        % Extract raw data
        X = Vdc(:,1,i);
        Y = dIdV_nor(:,1,i);

        % Bias voltage in mV (rescaled)
        xdata(:,j,i) = 1000 .* X;        
        ydata(:,j,i) = Y;

        % Define fitting range
        [~, Lindex] = min(abs(xdata - xRange(1)));
        [~, Rindex] = min(abs(xdata - xRange(2)));

        % Choose fitting model
        if isNodal
            [fitParam(:,j,i),resnorm,residual,exitflag,output,lambda,Jacob] = ...
                lsqcurvefit(@DynesFit_nodal_gauss_fraction_Base, ...
                param0, xdata(Lindex:Rindex,j,i), ydata(Lindex:Rindex,j,i), ...
                ParamLowerBound, ParamUpperBound, options);
        else
            [fitParam(:,j,i),resnorm,residual,exitflag,output,lambda,Jacob] = ...
                lsqcurvefit(@DynesFit_swave_gauss_fraction_Base, ...
                param0, xdata(Lindex:Rindex,j,i), ydata(Lindex:Rindex,j,i), ...
                ParamLowerBound, ParamUpperBound, options);
        end

        % Confidence intervals for fit parameters
        ParamErr(:,:,j,i) = nlparci(fitParam(:,j,i)',residual,'jacobian',Jacob);

        % Display extracted delta and gamma with errors
        parameters = {'delta (meV)'; 'gamma (meV)'};
        value = fitParam(1:2,j,i);
        min_err(:,j,i) = [(ParamErr(1,1,j,i)-fitParam(1,j,i)); ...
                          (ParamErr(2,1,j,i)-fitParam(2,j,i))];
        max_err(:,j,i) = [(ParamErr(1,2,j,i)-fitParam(1,j,i)); ...
                          (ParamErr(2,2,j,i)-fitParam(2,j,i))];
        disp(table(parameters,value,squeeze(min_err(:,j,i)),squeeze(max_err(:,j,i))))
    end
end

%% Plot: Measured Data with Fits
figure(1023); clf;
for i = 1:length(Vtglist)
    for j = 1:fitRange
        % Plot measured spectra
        plot(xdata(:,j,i), dIdV_nor(:,1,i) + 0.3*(j-1), ...
            'LineWidth', 1.5, 'Color', lineStyles(j,:), ...
            'DisplayName',['T = ' num2str(round(T_ave(j),2)) ' K']);
        hold on; xlim([-0.5, 0.5]);

        % Plot fitted curve
        if isNodal
            yydata(:,j,i) = DynesFit_nodal_gauss_fraction_Base(fitParam(:,j,i), xdata(:,j,i));
            title(['Nodal fitting with Gaussian at V_{tg} = ' num2str(Vtglist(i)) ' V'])
        else
            yydata(:,j,i) = DynesFit_swave_gauss_fraction_Base(fitParam(:,j,i), xdata(:,j,i));
            title(['s-wave fitting with Gaussian at V_{tg} = ' num2str(Vtglist(i)) ' V'])
        end
        plot(xdata(Lindex:Rindex,j,i), yydata(Lindex:Rindex,j,i) + 0.3*(j-1), ...
            '-.','LineWidth',1,'Color','black', ...
            'DisplayName',['Fit, T = ' num2str(round(T_ave(j),2)) ' K']);
    end
    xlabel('V_{dc} (mV)');
    ylabel('dI/dV_{norm} (a.u.)');
    ylim([0 1.4])
end

%% Plot: Gaussian Gap Distribution
figure(151); clf;
for i = 1:length(Vtglist)
    for j = 1:fitRange
        sigma = fitParam(5,j,i);
        delta0 = fitParam(1,j,i);
        delta = -2:0.001:2; % Gap axis (meV)

        % Gaussian distribution
        Gauss = @(d) 1./(sqrt(2*pi).*sigma) .* exp(-(d-delta0).^2/(2*sigma.^2));

        plot(delta, Gauss(delta), 'Color', lineStyles(j,:), ...
            'DisplayName',['T = ' num2str(round(T_ave(j),2)) ' K']);
        hold on
    end
    xlabel('Superconducting gap (meV)');
    ylabel('Probability density function');
    if isNodal
        title(['Nodal gap distribution at V_{tg} = ' num2str(Vtglist(i)) ' V']);
    else
        title(['s-wave gap distribution at V_{tg} = ' num2str(Vtglist(i)) ' V']);
    end
    legend('show')
end

%% Plot: Temperature Dependence of Fit Parameters
figure(152); clf;
for i = 1:length(Vtglist)
    plot(T_ave(1:fitRange), fitParam(1,:,i), '-o','Color','b','DisplayName','Gap');
    hold on
    plot(T_ave(1:fitRange), fitParam(2,:,i), '-o','Color','r','DisplayName','Gamma');
    plot(T_ave(1:fitRange), fitParam(5,:,i), '-o','Color','m','DisplayName','Sigma');

    xlabel('T (K)');
    ylabel('Fit parameters (meV)');
    if isNodal
        title('Nodal parameter temperature dependence');
    else
        title('s-wave parameter temperature dependence');
    end
    legend('show','Location','northwest')
end


